function [k_c, Ma_c] = neutral_stab_curves(k0)

N = 200;
N_search = 30;
ds = 1;

tol = 1e-8;

% preallocate
Ma_c = zeros(N, 1);
k_c = Ma_c;

ev_hat = zeros(N_search, 1);
k_hat = ev_hat;
Ma_hat = ev_hat;
theta = ev_hat;

p = params;


% first and second computation using secant method
for n = 1:2
    
    if n == 2
        p.Ma = 0.9 * Ma_c(1);
    end
    
    for its = 1:N_search
        
        if (its == 1)
            if (n == 1)
                k_hat(its) = k0;
            else
                k_hat = k_hat(n-1);
            end
        elseif (its == 2)
            k_hat(its) = k_hat(its-1) + 1e-5;
        else
            k_hat(its) = k_hat(its - 1) - ev_hat(its-1) / dev;
        end
        
        ev_hat(its) = comp_eigs(k_hat(its), p);
        
        if (abs(ev_hat(its)) < tol)
            k_c(n) = k_hat(its);
            Ma_c(n) = p.Ma;
            break
            
        else
            
            if (its > 1)
                dev = (ev_hat(its) - ev_hat(its-1)) / (k_hat(its) - k_hat(its-1));
            end
            
        end        
        
    end
end

conv = 1;
n = n + 1;
while (n <= N)
    
    if (conv == 1)
        ds = min(10 * ds, 1);
    else
        ds = ds / 10;
    end
    
    fprintf('n = %d\n', n);
    
    t_vec = [k_c(n-1) - k_c(n-2), Ma_c(n-1) - Ma_c(n-2)];
%     t_vec = t_vec / norm(t_vec);
    n_vec = [t_vec(2), -t_vec(1)] / norm(t_vec);
    
    k_hat(1) = k_c(n-1) + ds * t_vec(1);
    Ma_hat(1) = Ma_c(n-1) + ds * t_vec(2);
    
    p.Ma = Ma_hat(1);
    
    theta(1) = 0;
    ev_hat(1) = comp_eigs(k_hat(1), p);
    if abs(ev_hat(1)) < tol
        
        k_c(n) = k_hat(1);
        Ma_c(n) = Ma_hat(1);
        conv = 1;
        n = n + 1;
        
    else
        
        for its = 2:N_search
            
            if (its == 2)
                theta(2) = 1e-5;
            else
                theta(its) = theta(its-1) - ev_hat(its-1) / dev_dtheta;
            end
            
            k_hat(its) = k_hat(1) + theta(its) * n_vec(1);
            Ma_hat(its) = Ma_hat(1) + theta(its) * n_vec(2);
            
            p.Ma = Ma_hat(its);
            ev_hat(its) = comp_eigs(k_hat(its), p);
            
            fprintf('%.4e\n', abs(ev_hat(its)));
            
            if (abs(ev_hat(its)) < tol)
                k_c(n) = k_hat(its);
                Ma_c(n) = Ma_hat(its);
                fprintf('-----------------------\n');
                conv = 1;
                n = n + 1;
                break
            end
            
            dev_dtheta = (ev_hat(its) - ev_hat(its-1)) / (theta(its) - theta(its-1));
            
        end
        
        if its == N_search
            fprintf('no convergence!\n');
            conv = 0;
%             k_c = k_c(1:n-1);
%             Ma_c = Ma_c(1:n-1);
%             break
        end
        
    end
    
end


% plot(k_c, Ma_c, 'k*');


